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Finite-difference,  time-domain  (FDTD)  calculations  are  typically  performed  with  partial  differential 
equations  that  are  first  order  in  time.  Equation  sets  appropriate  for  FDTD  calculations  in  a  moving 
inhomogeneous  medium  (with  an  emphasis  on  the  atmosphere)  are  derived  and  discussed  in  this 
paper.  Two  candidate  equation  sets,  both  derived  from  linearized  equations  of  fluid  dynamics,  are 
proposed.  The  first,  which  contains  three  coupled  equations  for  the  sound  pressure,  vector  acoustic 
velocity,  and  acoustic  density,  is  obtained  without  any  approximations.  The  second,  which  contains 
two  coupled  equations  for  the  sound  pressure  and  vector  acoustic  velocity,  is  derived  by  ignoring 
terms  proportional  to  the  divergence  of  the  medium  velocity  and  the  gradient  of  the  ambient 
pressure.  It  is  shown  that  the  second  set  has  the  same  or  a  wider  range  of  applicability  than  equations 
for  the  sound  pressure  that  have  been  previously  used  for  analytical  and  numerical  studies  of  sound 
propagation  in  a  moving  atmosphere.  Practical  FDTD  implementation  of  the  second  set  of  equations 
is  discussed.  Results  show  good  agreement  with  theoretical  predictions  of  the  sound  pressure  due  to 
a  point  monochromatic  source  in  a  uniform,  high  Mach  number  flow  and  with  Fast  Field  Program 
calculations  of  sound  propagation  in  a  stratified  moving  atmosphere.  ©  2005  Acoustical  Society  of 
America.  [DOI;  10.1121/1.1841531] 

PACS  numbers;  43.20.Bi,  43.28. Js  [MO]  Pages;  503-517 


I.  INTRODUCTION 

Finite-difference,  time-domain  (FDTD)  techniques  have 
drawn  substantial  interest  recently  due  to  their  ability  to 
readily  handle  complicated  phenomena  in  outdoor  sound 
propagation  such  as  scattering  from  buildings  and  trees,  dy¬ 
namic  turbulence  fields,  complex  moving  source  distribu¬ 
tions,  and  propagation  of  transient  signals. These  phenom¬ 
ena  are  difficult  to  handle  with  frequency-domain 
techniques  that  are  currently  widely  used,  such  as  parabolic 
equation  approximations  and  the  Fast  Field  Program  (FFP). 
FDTD  techniques  typically  solve  coupled  sets  of  partial  dif¬ 
ferential  equations  that  are  first  order  in  time.  In  this  regard, 
they  are  a  departure  from  methodologies  such  as  the  para¬ 
bolic  approximation,  which  solve  a  single  equation  for  the 
sound  pressure  that  is  second  order  in  time.  Many  such 
single  equations  for  the  sound  pressure  in  a  moving  inhomo¬ 
geneous  medium  are  known  in  the  literature  (see  Refs.  9-14 


“'Portions  of  this  work  were  presented  in  V.  E.  Ostashev,  L.  Liu,  D.  K. 
Wilson,  M.  L.  Moran,  D.  F.  Aldridge,  and  D.  Marlin,  “Starting  equations 
for  direct  numerical  simulation  of  sound  propagation  in  the  atmosphere,” 
Proceedings  of  the  10th  International  Symposium  on  Long  Range  Sound 
Propagation,  Grenoble,  France,  Sept.  2002,  pp.  73-81. 


and  references  therein).  Although  these  equations  were  ob¬ 
tained  with  different  assumptions  and/or  approximations,  all 
contain  second-  or  higher-order  derivatives  of  the  sound 
pressure  with  respect  to  time,  and  are  therefore  not  amenable 
to  first-order  FDTD  techniques.  Our  main  goal  in  the  present 
paper  is  to  derive  equation  sets  that  are  appropriate  as  start¬ 
ing  equations  in  FDTD  simulations  of  sound  propagation  in  a 
moving  inhomogeneous  atmosphere  and  to  study  the  range 
of  applicability  of  these  sets. 

The  most  general  possible  approach  to  sound  propaga¬ 
tion  in  a  moving  inhomogeneous  medium  would  be  based  on 
a  direct  solution  of  the  complete  set  of  linearized  equations 
of  fluid  dynamics,®  "  which  are  first-order  partial  differ¬ 
ential  equations.  Although  this  set  could  be  used  as  starting 
equations  for  FDTD  codes,  even  with  modern  computers  it  is 
too  involved  to  be  practical.  Furthermore,  this  set  contains 
the  ambient  pressure  and  entropy,  which  are  not  usually  con¬ 
sidered  in  studies  of  sound  propagation  in  the  atmosphere. 
Therefore,  it  is  worthwhile  to  And  simplified  equation  sets 
for  use  in  FDTD  calculations. 

In  the  present  paper,  the  complete  set  of  linearized  equa¬ 
tions  of  fluid  dynamics  in  a  moving  inhomogeneous  medium 
is  reduced  to  two  simpler  sets  that  are  first  order  in  time  and 
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amenable  to  FDTD  implementation.  The  first  set  contains 
three  coupled  equations  involving  the  sound  pressure,  vector 
acoustic  (particle)  velocity,  and  acoustic  density.  No  approxi¬ 
mations  are  made  in  deriving  this  set.  The  second  set  con¬ 
tains  two  coupled  equations  for  the  sound  pressure  and  vec¬ 
tor  acoustic  velocity.  Although  the  second  set  describes 
sound  propagation  only  approximately,  the  assumptions  in¬ 
volved  in  deriving  the  second  set  are  quite  reasonable  in 
atmospheric  acoustics:  Terms  proportional  to  the  divergence 
of  the  medium  velocity  and  the  gradient  of  the  ambient  at¬ 
mospheric  pressure  are  ignored.  To  better  understand  the 
range  of  applicability  of  the  second  set,  we  compare  the  set 
with  equations  for  the  sound  pressure  that  have  been  previ¬ 
ously  used  in  analytical  and  numerical  studies  of  sound 
propagation  in  a  moving  atmosphere.  It  is  shown  that  the 
second  set  has  the  same  or  a  wider  range  of  applicability  than 
these  equations  for  the  sound  pressure. 

Furthermore  in  the  present  paper,  a  basic  numerical  al¬ 
gorithm  for  solving  the  second  set  of  equations  in  two- 
dimensional  (2-D)  moving  inhomogeneous  media  is  devel¬ 
oped.  Issues  related  to  the  finite-difference  approximation  of 
the  spatial  and  temporal  derivatives  are  discussed.  FDTD 
solutions  are  obtained  for  a  homogenous  uniformly  moving 
medium  and  for  a  stratified  moving  atmosphere.  The  first  of 
these  solutions  is  compared  with  an  analytical  formula  for 
the  sound  pressure  due  to  a  point  monochromatic  source  in  a 
uniformly  moving  medium.  The  second  solution  is  compared 
with  predictions  from  before  FFR 

Although  the  explicit  emphasis  of  the  discussion  in  this 
paper  is  on  sound  propagation  in  a  moving  inhomogeneous 
atmosphere,  most  of  the  derived  equations  are  also  valid  for 
a  general  case  of  sound  propagation  in  a  moving  inhomoge¬ 
neous  medium  with  an  arbitrary  equation  of  state,  e.g.,  in  the 
ocean  with  currents.  Equations  presented  in  the  paper  are 
also  compared  with  those  known  in  aeroacoustics. 

The  paper  is  organized  as  follows.  In  Sec.  II,  we  con¬ 
sider  the  complete  set  of  equations  of  fluid  dynamics  and 
their  linearization.  In  Sec.  Ill,  the  linearized  equations  are 
reduced  to  the  set  of  three  coupled  equations  for  the  sound 
pressure,  acoustic  velocity,  and  acoustic  density.  In  Sec.  IV, 
we  consider  the  set  of  two  coupled  equations  for  the  acoustic 
pressure  and  acoustic  velocity.  Numerical  implementation  of 
this  set  is  considered  in  Sec.  V. 


II.  EQUATIONS  OF  FLUID  DYNAMICS  AND  THEIR 
LINEARIZATION 

Let  P(R,f)  be  the  pressure,  E(RT)  the  density,  "^(R,?) 
the  velocity  vector,  and  S(R,t)  the  entropy  in  a  medium. 
Here,  R=(x,y,z)  are  the  Cartesian  coordinates,  and  t  is 
time.  These  functions  satisfy  a  complete  set  of  fluid  dynamic 
equations  (e.g.  Ref.  16): 


d  \  VP 

— -FvV  v+  — -g=F/e, 

1  Q 

(1) 

d  \  _ 

—  +  yV\Q  +  QV-y=QQ, 

(2) 

p=pie,s).  (4) 

In  Eqs.  (l)-(4),  V  =  {dldx,dldy,d/dz),  g=(0,0,g)  is  the  ac¬ 
celeration  due  to  gravity,  and  F  and  Q  characterize  a  force 
acting  on  the  medium  and  a  mass  source,  respectively.  Eor 
simplicity,  we  do  not  consider  the  case  when  a  passive  com¬ 
ponent  is  dissolved  in  a  medium  (e.g.,  water  vapor  in  the  dry 
air,  or  salt  in  water).  This  case  is  considered  in  detail 
elsewhere.®’*^ 

If  a  sound  wave  propagates  in  a  medium,  in  Eqs.  (l)-(4) 
P,  Q,  V,  and  S  can  be  expressed  in  the  following  form:  P 
=  P+p,  Q=Q+  rj,  v  =  v-l-w,  and  S  =  S  +  s.  Here,  P,  Q,  v, 
and  S  are  the  ambient  values  (i.e.,  the  values  in  the  absence 
of  a  sound  wave)  of  the  pressure,  density,  medium  velocity, 
and  entropy  in  a  medium,  and  p,  77,  w,  and  s  are  their  fluc¬ 
tuations  due  to  a  propagating  sound  wave.  In  order  to  obtain 
equations  for  a  sound  wave,  Eqs.  (l)-(4)  are  linearized  with 
respect  to  p,  77,  w,  and  s.  Assuming  that  a  sound  wave  is 
generated  by  the  mass  source  Q  and/or  the  force  F  and  in¬ 
troducing  the  full  derivative  with  respect  to  time  didt 
=  d! dt  +  y'^ ,  we  have 


c/w  Vp  7]V  P 

(5) 

-F(wV)v-F  2 

dt  ^  Q  E 

d  77 

-F  (wV)g  -F  g  V-w-F  rjV ■\=  QQ, 

(6) 

ds 

— -F(wV)5  =  0, 
dt  ^  ’ 

(7) 

p=  rjc^  +  hs. 

(8) 

Here,  c=  \ldP(Q,S)/dg  is  the  adiabatic  sound  speed,  and 
the  parameter  h  is  given  hy  h  =  dP{Q,S)/dS.  The  set  of  Eqs. 
(5)-(8)  provides  a  most  general  description  of  sound  propa¬ 
gation  in  a  moving  inhomogeneous  medium  with  only  one 
component.  In  order  to  calculate  p,  77,  w,  and  s,  one  needs  to 
know  the  ambient  quantities  c,  Q,  v,  P,  S,  and  h.  Note  that 
Eqs.  (5)-(8)  describe  the  propagation  of  both  acoustic  and 
internal  gravity  waves,  as  well  as  vorticity  and  entropy 
waves  (e.g..  Ref.  18). 

Equations  (5)-(8)  were  derived  for  the  first  time  by 
Blokhintzev  in  1946.^^  Since  then,  these  equations  have  been 
widely  used  in  studies  of  sound  propagation  (e.g..  Refs. 
9-11).  In  the  general  case  of  a  moving  inhomogeneous  me¬ 
dium,  Eqs.  (5)-(8)  cannot  be  exactly  reduced  to  a  single 
equation  for  the  sound  pressure  p.  In  the  literature,  Eqs. 
(5)-(8)  have  been  reduced  to  equations  for  p,  making  use  of 
different  approximations  or  assumptions  about  the  ambient 
medium.  These  equations  for  p  were  subsequently  used  for 
analytical  and  numerical  studies  of  sound  propagation.  They 
are  discussed  in  Sec.  IV.  Note  that  the  equations  for  p  known 
in  the  literature  contain  the  following  ambient  quantities:  c, 
Q,  and  V.  On  the  other  hand,  the  linearized  equations  of  fluid 
dynamics,  Eqs.  (5)-(8),  contain  not  only  c,  Q,  and  v,  but 
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also  P,  S,  and  h.  This  fact  indicates  that  the  effect  of  P,  S, 
and  h  on  sound  propagation  is  probably  small  for  most  of 
problems  considered  so  far  in  the  literature. 

The  effect  of  medium  motion  on  sound  propagation  is 
also  studied  in  aeroacoustics,  e.g.,  see  Refs.  12,  19-24  and 
references  therein.  In  aeroacoustics,  the  starting  equations 
coincide  with  Eqs.  (l)-(4)  but  might  also  include  terms  de¬ 
scribing  viscosity  and  thermal  conductivity  in  a  medium.  Us¬ 
ing  these  equations  of  fluid  dynamics,  equations  for  sound 
waves  are  derived  which  have  some  similarities  with  Eqs. 
(5) -(8).  Eor  example,  Eqs.  (5) -(8)  are  equivalent  to  Eqs. 
(1.11)  from  Ref.  12,  and  Eq.  (6)  can  be  found  in  Refs.  19,  22, 
23.  The  main  difference  between  Eqs.  (5)-(8)  and  those  in 
aeroacoustics  are  sound  sources.  In  atmospheric  acoustics,  in 
Eqs.  (5)-(8)  the  sources  F  and  Q  are  assumed  to  be  known 
and  are  loudspeakers,  car  engines,  etc.  In  aeroacoustics,  these 
sources  have  to  be  calculated  and  are  those  due  to  ambient 
flow.  Eurthermore  in  some  formulations  in  aeroacoustics,  the 
left-hand  side  of  Eq.  (5)  contains  nonlinear  terms. Note 
that  EDTD  calculations  are  nowadays  widely  used  in  aeroa¬ 
coustics,  e.g..  Refs.  19,  20,  24. 

Also  note  that  in  aeroacoustics  it  is  sometimes  assumed 
that  the  ambient  medium  is  incompressible  and/or  isentropic, 
i.e.,  5  =  const.  Generally,  these  assumptions  are  inappropriate 
for  atmospheric  acoustics.  Indeed,  sound  waves  can  be  sig¬ 
nificantly  scattered  by  density  fluctuations,  e.g.,  see  Sec. 
6.1.4  from  Ref.  9.  Eurthermore,  in  a  stratified  atmosphere  S 
depends  on  the  height  above  the  ground.  The  range  of  appli¬ 
cability  of  the  assumption  5  =  const  (which  is  equivalent  to 
s  =  0  ox  p  =  c^rf)  is  studied  in  Sec.  2.2.4  from  Ref.  9.  Eor  a 
stratified  medium,  this  assumption  is  not  applicable  if  the 
scale  of  the  ambient  density  variations  is  smaller  than  the 
sound  wave  length  or  if  the  ambient  density  noticeably 
changes  with  height. 


III.  SET  OF  THREE  COUPLED  EQUATIONS 

A.  Moving  medium  with  an  arbitrary  equation  of  state 

Applying  the  operator  {d! dt  +  y-'^)  to  both  sides  of  Eq. 
(4)  and  using  Eq.  (3),  we  have 

l+,.v)p  =  c2(^  +  v.v)g,  (9) 

where  c^=  dP{Q,S)ldQ  differs  from  the  square  of  the  adia¬ 
batic  sound  speed  c^  =  dP(Q,S)ldQ.  Using  Eq.  (2),  Eq.  (9) 
can  be  written  as 


— -l-vV  \P  +  c^QV-y=c^QQ. 
at 


(10) 


The  next  step  is  to  linearize  Eq.  (10)  to  obtain  an  equation 
for  acoustic  quantities.  To  do  so  we  need  to  calculate  the 
value  of  c^=  dP(g,S)/dQ  to  the  first  order  in  acoustic  per¬ 
turbations.  In  this  formula,  we  express  Q  and  S  as  the  sums 
Q  =  g+  and  S  =  S  +  s,  decompose  the  function  P  into  Tay¬ 
lor  series,  and  keep  the  terms  of  the  first  order  in  77  and  s: 


dP{e,s) 

dQ 


—  PiQ+rf,S  +  s) 
dQ 


d 

dP{Q,S) 


dP{Q,S)  dP{Q,S) 
P{e,S)P - - r]+ - - s 


dQ 


dS 


dQ 


d^P(Q,S)  d^P{Q,S) 

3 - — - ’7-1 - — S. 


dQ^ 


dQdS 


(11) 


The  first  term  in  the  last  line  of  this  equation  is  equal  to  c^. 
Denoting  /3=d^P{Q,S)ldQ^  and  a=  d^P{Q,S)/dQdS,  we 
have  c^  =  c^  +  PyP  as  =  c^  +  {c^)' .  Here,  (c^)'  =  ^7/+  as 
are  fluctuations  in  the  squared  sound  speed  due  to  a  propa¬ 
gating  sound  wave.  In  this  formula,  s  can  be  replaced  by  its 
value  from  Eq.  (8);  s  =  (p  —  c^rf)/h.  As  a  result,  we  obtain 
the  desired  formula  for  fluctuations  in  the  squared  sound 
speed;  (c^)'  =  {p—  ac^lh)7j+  aplh. 

Now  we  can  linearize  Eq.  (10).  In  this  equation,  we 
express  P,  Q,  v,  and  as  the  sums;  P  =  P+p,  Q  =  Q+  p, 
v=v-l-w,  and  c^  =  c^  +  (c^)'.  Linearizing  the  resulting  equa¬ 
tion  with  respect  to  acoustic  quantities,  we  have 


—  -f  gc^V-w-l- wVf’-l-(c^77-l-  g(c^)')V-v=  Qc^Q- 

(12) 

In  this  equation,  (c^)'  is  replaced  by  its  value  obtained 
above.  As  a  result,  we  arrive  at  the  following  equation  for 
dpidt: 


dp 

dt 


+  gc^V-w-l- wVP-l-{[g;S-l-c^(  1  —  ag//?)]  77 


+  {aQlh)p}^  ■y=  Qc^Q.  (13) 

Equations  (5),  (6),  and  (13)  comprise  a  desired  set  of 
three  coupled  equations  for  p,  w,  and  77.  This  set  was  ob¬ 
tained  from  linearized  equations  of  fluid  dynamics,  Eqs.  (5)- 
(8),  without  any  approximations.  The  set  can  be  used  as  start¬ 
ing  equations  for  EDTD  simulations.  In  this  set,  one  needs  to 
know  the  following  ambient  quantities;  c,  Q,  v,  P,  a,  j3,  and 
h. 


B.  Set  of  three  equations  for  an  ideal  gas 

In  most  applications,  the  atmosphere  can  be  considered 
as  an  ideal  gas.  In  this  case,  the  equation  of  state  reads  (e.g.. 
Refs.  9,  17)  as 

^  =  ^o(e/eo)"exp[(7-l)(5-5o)/Rj,  (14) 

where  y=  1.4  is  the  ratio  of  specific  heats  at  constant  pres¬ 
sure  and  constant  volume,  is  the  gas  constant  for  the  air, 
and  the  subscript  0  indicates  reference  values  of  P,  Q,  and  S. 
Using  Eq.  (14),  the  sound  speed  c  and  the  coefficients  a,  j3, 
and  h  appearing  in  Eq.  (13)  can  be  calculated;  c^=yP/Q, 
a=  y(y— l)P/(QRg),  /3=y(y—l)PlQ^,  and  h  =  (y 
—  l)PIRa  ■  Substituting  these  values  into  Eq.  (13),  we  have 

+  Qc^^  ■yx  +  yx-'^  P  +  yp'^  ■y=  Qc^Q.  (15) 

A  set  of  Eqs.  (5),  (6),  and  (15)  is  a  closed  set  of  three 
coupled  equations  for  p,  w,  and  77  for  the  case  of  an  ideal 
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gas.  To  solve  these  equations,  one  needs  to  know  the  follow¬ 
ing  ambient  quantities:  c,  Q,  v,  and  P. 

Let  us  compare  Eqs.  (5),  (6),  and  (15)  with  a  closed  set 
of  equations  for  p  and  w  from  Ref.  1;  see  Eqs.  (12)  and  (13) 
from  that  reference.  The  latter  set  was  used  in  Refs.  1,  2  as 
starting  equations  for  EDTD  simulations  of  outdoor  sound 
propagation.  If  Q  =  0,  Eq.  (15)  in  the  present  paper  is  essen¬ 
tially  the  same  as  Eq.  (13)  from  Ref.  1.  [Note  that  Eq.  (15)  is 
also  used  in  aeroacoustics,  e.g.,  Ref.  19.]  Eurthermore  for  the 
case  of  a  nonabsorbing  medium,  Eq.  (12)  from  Ref.  1  is 
given  by 

dw  Vp  pVP 

— -(-(wV)v-l - ^=0-  (16) 

dt  Q  yPQ 

Let  us  show  that  this  equation  is  an  approximate  version  of 
Eq.  (5)  in  the  present  paper.  Indeed,  in  Eq.  (5)  we  replace  rj 
by  its  value  from  Eq.  (8):  rj  =  (p  —  hs)lc^,  and  assume  that 
s  =  0.  If  F  =  0,  the  resulting  equation  coincides  with  Eq.  (16). 
Thus,  for  an  ideal  gas  and  F  =  0  and  Q  =  Q,  Eqs.  (12)  and 
(13)  from  Ref.  1  are  equivalent  to  Eqs.  (5)  and  (15)  in  the 
present  paper  if  s  can  be  set  to  0.  The  range  of  applicability 
of  the  approximation  s  =  0  is  considered  above. 

IV.  SET  OF  TWO  COUPLED  EQUATIONS 
A.  Set  of  equations  for  p  and  w 

In  atmospheric  acoustics,  Eqs.  (5)  and  (13)  can  be  sim- 
plihed  since  v  is  always  much  less  than  c.  Eirst,  using  Ref. 
16,  it  can  be  shown  that  V-v~v^l{c^L),  where  L  is  the 
length  scale  of  variations  in  the  density  Q.  Therefore,  in  Eq. 
(13)  the  term  proportional  to  V-v  can  be  ignored  to  order 
v^lc^ ■  Second,  in  Eqs.  (5)  and  (13)  the  terms  proportional  to 
VP  can  also  be  ignored.  Indeed,  in  a  moving  inhomoge¬ 
neous  atmosphere  VP  is  of  the  order  v^lc^  so  that  these 
terms  can  be  ignored  to  order  vie.  Eurthermore,  in  a  strati- 
hed  atmosphere,  ^P=—gQ,  where  g  is  the  acceleration  due 
to  gravity.  It  is  known  that,  in  linearized  equations  of  fluid 
dynamics,  terms  proportional  to  g  are  important  for  internal 
gravity  waves  and  can  be  omitted  for  acoustic  waves. 

With  these  approximations,  Eqs.  (13)  and  (5)  become 


-Fvv|/7-Fec^V-w=ec^2, 

(17) 

\  Vp 

-FvVj  w-F(wV)v-F  —  =  F/e. 

(18) 

Equations  (17)  and  (18)  comprise  the  desired  closed  set  of 
two  coupled  equations  for  p  and  w.  This  set  can  also  be  used 
in  EDTD  simulations  of  sound  propagation  in  the  atmo¬ 
sphere.  In  order  to  solve  this  set,  one  needs  to  know  the 
following  ambient  quantities:  c,  Q,  and  v.  These  ambient 
quantities  appear  in  equations  for  the  sound  pressure  p  that 
have  been  most  often  used  for  analytical  and  numerical  stud¬ 
ies  of  sound  propagation  in  moving  media.  The  set  of  Eqs. 
(17)  and  (18)  is  simpler  than  the  set  of  three  coupled  equa¬ 
tions,  Eqs.  (5),  (6),  and  (13),  and  does  not  contain  the  ambi¬ 
ent  quantities  P,  a,  (3,  and  h.  It  can  be  shown  that  Eqs.  (17) 
and  (18)  describe  the  propagation  of  acoustic  and  vorticity 
waves  but  do  not  describe  entropy  or  internal  gravity  waves. 


Equations  (17)  and  (18)  were  derived  in  Ref.  25  [see 
also  Eqs.  (2.68)  and  (2.69)  from  Ref.  9]  using  a  different 
approach.  In  these  references,  Eqs.  (17)  and  (18)  were  de¬ 
rived  for  the  case  of  a  moving  inhomogeneous  medium  with 
more  than  one  component  (e.g.,  humid  air  or  salt  water). 
Equations  (17)  and  (18)  are  somewhat  similar  to  the  starting 
equations  in  EDTD  simulations  used  in  Ref.  3;  see  Eqs.  (10) 
and  (12)  from  that  reference.  The  last  of  these  equations 
coincides  with  Eq.  (17)  while  the  first  is  given  by 

dVi  Vp 

- wX(VXv)-l - l-V[wv]  =  0.  (19) 

dt  Qq 

Using  vector  algebra,  the  left-hand  side  of  this  equation  can 
be  written  as  a  left-hand  side  of  Eq.  (18)  plus  an  extra  term 
vX(VXw).  Equations  (10)  and  (12)  from  Ref.  3  were  ob¬ 
tained  using  several  assumptions  that  were  not  employed  in 
the  present  paper  when  deriving  Eqs.  (17)  and  (18):  d\ldt 
=  dgl dt  =  V Q  =  0,  c  is  constant,  and  dwl dt>\X{V Xw) .  It 
follows  from  the  last  inequality  that  the  “extra”  term 
vX(VXw)  in  Eq.  (19)  can  actually  be  omitted.  Note  that  in 
Ref.  4  different  starting  equations  were  used  in  simulations 
of  sound  propagation  in  a  muffler  with  a  low  Mach  number 
flow.  The  use  of  Eq.  (19)  resulted  in  increase  of  stability  in 
such  simulations. 

Also  note  that  equations  for  p  and  w  similar  to  Eqs.  (17) 
and  (18)  are  used  in  aeroacoustics,  e.g.  Refs.  20,  24.  The 
left-hand  sides  of  Eqs.  (7)  in  Ref.  20  contain  several  extra 
terms  in  comparison  with  the  left-hand  sides  of  Eqs.  (17)  and 
(18)  which,  however,  vanish  if  VP  =  0  and  V-v=0.  The  left- 
hand  sides  of  Eqs.  (75)  and  (76)  in  Ref.  24  also  contain  extra 
terms  in  comparison  with  the  left-hand  sides  of  Eqs.  (17)  and 
(18),  e.g.,  terms  proportional  to  the  gradients  of  c  and  Q.  The 
right-hand  sides  of  the  equations  in  Refs.  20,  24  describe 
aeroacoustic  sources  and  differ  from  those  in  Eqs.  (17)  and 
(18). 

At  the  beginning  of  this  section,  we  provided  sufficient 
conditions  for  the  applicability  of  Eqs.  (17)  and  (18).  Actu¬ 
ally,  the  range  of  applicability  of  these  equations  can  be 
much  wider.  Note  that  it  is  quite  difficult  to  estimate  with 
what  accuracy  one  can  ignore  certain  terms  in  differential 
equations.  We  will  study  the  range  of  applicability  of  Eqs. 
(17)  and  (18)  by  comparing  them  with  equations  for  the 
sound  pressure  p  presented  in  Secs.  IVB-IVE,  which  have 
been  most  often  used  for  analytical  and  numerical  studies  of 
sound  propagation  in  moving  media  and  whose  ranges  of 
applicability  are  well  known.  This  will  allow  us  to  show  that 
Eqs.  (17)  and  (18)  have  the  same  of  a  wider  range  of  appli¬ 
cability  than  these  equations  for  p  and,  in  many  cases,  de¬ 
scribe  sound  propagation  to  any  order  inv/c.  Eor  simplicity, 
in  the  rest  of  this  section,  we  assume  that  F  =  0,  2  =  0,  and 
the  medium  velocity  is  subsonic. 

B.  Nonmoving  medium 

Consider  the  case  of  a  nonmoving  medium  when  v=0. 
In  this  case,  the  set  of  linearized  equations  of  fluid  dynamics, 
Eqs.  (5)-(8),  can  be  exactly  (without  any  approximations) 
reduced  to  a  single  equation  for  sound  pressure  p  (e.g.,  see 
Eq.  (1.11)  from  Ref.  11): 
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d 

lit 


1  dp\ 

Qc^  dt  I 


-V- 


=  0. 


(20) 


For  the  considered  case  of  a  nonmoving  medium,  Eqs. 
(17)  and  (18)  can  also  be  reduced  to  a  single  equation  for  p. 
This  equation  coincides  with  Eq.  (20).  Therefore,  Eqs.  (17) 
and  (18)  describe  sound  propagation  exactly  if  v=0. 


C.  Homogeneous  uniformly  moving  medium 


Wj^(r,z,f)  = 


dco  exp(/aT— i(uf)wj^(a,z,w). 

(27) 


Here,  r=  (x,y)  are  the  horizontal  coordinates,  a  is  the  hori¬ 
zontal  component  of  the  wave  vector,  co  is  the  frequency  of  a 
sound  wave,  and p,  w^,  and  are  the  spectral  densities  of 
p,  Wj. ,  and  .  We  substitute  Eqs.  (25)-(27)  into  Eqs.  (22)- 
(24).  As  a  result,  we  obtain  a  set  of  equations  for  p,  and 

Wj_: 


A  medium  is  homogeneous  and  uniformly  moving  if  the 
ambient  quantities  c,  v,  etc.  do  not  depend  on  R  and  t.  Eor 
such  a  medium,  the  linearized  equations  of  fluid  dynamics, 
Eqs.  (5)-(8)  can  also  be  exactly  reduced  to  a  single  equation 
for  p  (see  Sec.  2.3.6  from  Ref.  9  and  references  therein); 

(^  +  vv)  p-c^V^p  =  0.  (21) 

Eor  the  case  of  a  homogeneous  uniformly  moving  me¬ 
dium,  Eqs.  (17)  and  (18)  can  be  reduced  to  the  equation  for 
p  that  coincides  with  Eq.  (21).  Therefore,  Eqs.  (17)  and  (18) 
describe  sound  propagation  exactly  in  a  homogeneous  uni¬ 
formly  moving  medium.  In  particular,  they  correctly  account 
for  terms  of  any  order  in  v/c. 


D.  Stratified  moving  medium 

Now  let  us  consider  the  case  of  a  stratified  medium 
when  the  ambient  quantities  c,  Q,  etc.  depend  only  on  the 
vertical  coordinate  z-  We  will  assume  that  the  vertical  com¬ 
ponent  of  V  is  zero;  v=(vj^,0),  where  is  a  horizontal 
component  of  the  medium  velocity  vector.  In  this  subsection, 
we  reduce  Eqs.  (17)  and  (18)  to  a  single  equation  for  the 
spectral  density  of  the  sound  pressure  and  show  that  this 
equation  coincides  with  the  equation  for  the  spectral  density 
that  can  be  derived  from  Eqs.  (5)-(8). 

Eor  a  stratified  moving  medium,  Eq.  (17)  can  be  written 
as 


— Vl \p+Qc^\  V 


dw. 


— 1.0. 


(22) 


Here,  ^^  =  {8! dx,dl dy),  and  Wj^  and  are  the  horizontal 
and  vertical  components  of  the  vector  w=(w^  Equa¬ 
tion  (18)  can  be  written  as  two  equations; 


57+''^''^^ 


1  dp 

w^+  —  ^-  =  0, 

Q  dz 

(23) 

wx+H’,v;+  — =0- 

(24) 

Here,  =  Idz-  Let  p,  ,  and  be  expressed  as  Eou- 
rier  integrals; 


p{r,z,t) 


dw  exp{ia-r—  iwt)p{a,z,w). 


(25) 


w,{r,z,t) 


dw  exp{ia-r— iwt)w^{a,z,w), 

(26) 


■i{w  —  a-Vj^  )p  +  /  gc^a-Wj^  +  Qc^  =  0, 

8z 


(28) 


1  dp 

—  ;(&)  — a-Vi  )wy-\ - —  =  0, 

^  Q  8z 


(29) 


lap 

—  i{w  —  a-Vj^ )  -t-  v|vi>j.  H - =  0. 


(30) 


After  some  algebra,  this  set  of  equations  can  be  reduced  to  a 
single  equation  for  p: 


d-p  I  2a-v| 
8z^  \w-a-Vj^ 


g  I  dz 


{w-a-y^f 

- 9 - 0- 


P  =  0, 

(31) 


where  q'  =  dQldz. 

Eor  the  considered  case  of  a  stratified  moving  medium,  a 
single  equation  for  p  can  also  be  derived  from  Eqs.  (5)-(8) 
without  any  approximations.  This  equation  for  p  is  given  by 
Eq.  (2.61)  from  Ref.  9.  Setting  g  =  0  in  this  equation  (i.e. 
ignoring  internal  gravity  waves)  one  obtains  Eq.  (31).  There¬ 
fore,  Eqs.  (17)  and  (18)  describe  sound  propagation  exactly 
in  a  stratified  moving  medium,  and,  hence,  correctly  account 
for  terms  of  any  order  in  v/c. 


E.  Turbulent  medium 


Probably  the  most  general  of  the  equations  describing 
the  propagation  of  a  monochromatic  sound  wave  in  turbulent 
media  with  temperature  and  velocity  fluctuations  is  given  by 
Eq.  (6.1)  from  Ref.  9; 


A  -f  ^q(  1  -f  e)  ■ 


_  ,  2i  dVi  d^ 

Vln—  -V - ^ - 

Qq!  w  dXj  dXjdXj 


2ikn 

+ - ^vV 

Co 


p(R)  =  0. 


(32) 


Here,  A  =  d^/ dx^  + d^/ dy^  + d^/ dz^',  8  =  Co/c^— 1;  k^,  cq, 
and  go  ^6  the  reference  values  of  the  wave  number,  adia¬ 
batic  sound  speed,  and  density;  xi ,  X2,  X3  stand  for  x,  y,  z; 
Vi=Vjc,  V2  =  Vy,  V2  =  v^  are  the  components  of  the  medium 
velocity  vector  v;  and  repeated  subscripts  are  summed  from 
1  to  3.  Eurthermore,  the  dependence  of  the  sound  pressure  on 
the  time  factor  exp{  —  iwt)  is  omitted. 

The  range  of  applicability  of  Eq.  (32)  is  considered  in 
detail  in  Sec.  2.3  from  Ref.  9.  This  equation  was  used  for 
calculations  of  the  sound  scattering  cross  section  per  unit 
volume  of  a  sound  wave  propagating  in  a  turbulent  medium 
with  temperature  and  velocity  fluctuations.  Also  it  was  em¬ 
ployed  as  a  starting  equation  for  developing  a  theory  of  mul- 


J.  Acoust.  Soc.  Am.,  Vol.  117,  No.  2,  February  2005 


Ostashev  et  at.-.  Moving  media  finite  difference  time  domain  equations  507 


tiple  scattering  of  a  sound  wave  propagating  in  such  a  turbu¬ 
lent  medium;  see  Ref.  9  and  references  therein.  Furthermore, 
starting  from  Eq.  (32),  parabolic  and  wide-angle  parabolic 
equations  were  derived  and  used  in  analytical  and  numerical 
studies  of  sound  propagation  in  a  turbulent  medium,  e.g.. 
Ref  26.  For  example,  a  parabolic  equation  deduced  from  Eq. 
(32)  reads  as 


,  d0 

ecVi-V@-fpi-^  =  0, 


(40) 


d& 


v@ 


Wi  —;—+Pi - =  0. 

'  dt  Q 


(41) 


Equating  terms  proportional  to  k^,  we  obtain  another  set; 


dp 

2ikQ  —  +  ^^p  +  2kl 


1-f 


p  =  0. 


(33) 


Qc'^y/2-'^@+P2 


d& 

dt 


dp  I 
dt 


Qc^V-Wi, 


(42) 


Here,  the  predominant  direction  of  sound  propagation  coin¬ 
cides  with  the  x-axis,  =  dy'^,d~/ dz^),  and  £,^(,^,=  6 

-2vjca. 

In  Ref.  9,  Eq.  (32)  was  derived  starting  from  the  set  of 
Eqs.  (17)  and  (18)  and  using  some  approximations.  There¬ 
fore,  this  set  has  the  same  or  a  wider  range  of  applicability 
than  equations  for  p  that  have  been  used  in  the  literature  for 
analytical  and  numerical  studies  of  sound  propagation  in  a 
turbulent  medium  with  temperature  and  velocity  fluctuations. 


F.  Geometrical  acoustics 

Sound  propagation  in  a  moving  inhomogeneous  medium 
is  often  described  in  geometrical  acoustics  approximation 
which  is  applicable  if  the  sound  wavelength  is  much  smaller 
than  the  scale  of  medium  inhomogeneities.  In  geometrical 
acoustics,  the  phase  of  a  sound  wave  can  be  obtained  as  a 
solution  of  the  eikonal  equation,  and  its  amplitude  from  the 
transport  equation.  In  this  subsection,  starting  from  Eqs.  (17) 
and  (18),  we  derive  eikonal  and  transport  equations  and 
show  that  they  are  in  agreement  with  those  deduced  from 
Eqs.  (5)-(8). 

Let  us  express  p  and  w  in  the  following  form; 
p(R,f)  =  exp(/ko0(R,f))pA(RT),  (34) 

w(R,f)  =  exp(iko@(R,f))w^(R,f).  (35) 


d@  P2^0  d-^x  Vpi 

^  dt  Q  dt^'  Q 

From  Eq.  (41),  we  have 

Pi  V@ 

Wi  = -  (44) 

'  Q  d&/dt  ^  ^ 

Substituting  this  value  of  Wi  into  Eq.  (40),  we  obtain 

Pi  =  0.  (45) 


d& 


-I  -c^(V@)^ 


From  this  equation,  we  obtain  an  eikonal  equation  for  the 
phase  function; 


d& 

—  =  -c  V@  . 
dt 


(46) 


Here,  a  sign  in  front  of  |V0|  is  chosen  in  accordance  with  the 
time  convention  exp(  — ;wf).  Equation  (46)  coincides  exactly 
with  the  eikonal  equation  for  sound  waves  in  a  moving  in¬ 
homogeneous  medium  (e.g.,  see  Eq.  (3.15)  from  Ref.  9) 
which  can  be  derived  from  Eqs.  (5)-(8)  in  a  geometrical 
acoustics  approximation.  Thus,  in  this  approximation,  Eqs. 
(17)  and  (18)  exactly  describe  the  phase  of  a  sound  wave 
and,  hence,  account  for  terms  of  any  order  in  vie. 

Substituting  the  value  of  d&ldt  from  Eq.  (46)  into  Eq. 
(44),  we  have 


Here,  0(R,f)  is  the  phase  function,  and  p^  and  are  the 
amplitudes  of  p  and  w.  Substituting  Eqs.  (34)  and  (35)  into 
Eqs.  (17)  and  (18),  we  have 

/  T  d&\  dpjx 

ikol^  Qc^Wa  ■  V  0  +pa  -^j  =  -  V  •w,4 ,  (36) 


/  d% 


FxV@| 


dv/A 

dt 


-(w^- V)v- 


Q 

(37) 


In  geometrical  acoustics,  Pa  and  Wa  are  expressed  as  a  series 
in  a  small  parameter  proportional  to  l/kg ; 


Pa=Pi+j^  + 

IKq 


(38) 


Wo 

WA  =  Wi+Tf^  + 


(39) 


Substituting  Eqs.  (38)  and  (39)  into  Eqs.  (36)  and  (37)  and 
equating  terms  proportional  to  ko ,  we  arrive  at  a  set  of  equa¬ 
tions; 


Pi  V0  _pin 
Q  c|V0|  Qc  ’ 


(47) 


where  n=  V0/|V0|  is  the  unit  vector  normal  to  the  phase 
front.  Now  we  multiply  Eq.  (42)  by  d&ldt  and  multiply  Eq. 
(43)  by  c2gV0.  Then,  we  subtract  the  latter  equation  from 
the  former.  After  some  algebra  and  using  Eq.  (46),  it  can  be 
shown  that  the  sum  of  all  terms  proportional  to  p2  and  W2  is 
zero.  The  resulting  equation  reads  as 


dpx  dwi 

4-cn-Vpi -f  gen-  -f  gen- (wi  •  V)v 


-f  ge^V-Wi  =  0. 


(48) 


In  this  equation,  Wi  is  replaced  by  its  value  given  by  Eq. 
(47).  As  a  result,  we  obtain 

gn  d  lnpx\  1  dpx  n-Vpi  /npi\ 

- —  -  H - 9  — - 1 - f  g  V  ■  - 

c  dt\  Qc  j  c~  dt  c  \  Qc  j 


^  Pin-(n-V)v 

“t"  9 

c 


(49) 
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In  geometrical  acoustics,  the  amplitude  pj^  of  the  sound  pres¬ 
sure  is  approximated  by  p  j .  Equation  (49)  is  a  closed  equa¬ 
tion  for  Pi ;  i.e.,  it  is  a  transport  equation. 

The  second  term  on  the  left-hand  side  of  Eq.  (49)  can  be 
written  as 


1  dpi 

dt 


d  I pi\  Pi  dc^ 
dt\c^]  dt 


dt\c'^] 


,  Pi  PdQ 
dt 


(50) 


Here,  we  used  the  formula  dc^l dt  =  fSdQl dt\  see  Eq.  (2.63) 
from  Ref.  9.  According  to  Eq.  (2),  dQ/dt  in  Eq.  (50)  can  be 
replaced  with  —  gV-v.  When  deriving  Eqs.  (17)  and  (18), 
terms  proportional  to  V-v  were  ignored.  Therefore,  the  last 
term  on  the  right-hand  side  of  Eq.  (50)  should  also  be  ig¬ 
nored.  In  this  case,  Eq.  (49)  can  be  written  as 


gn  d  I  npj\  d 
c  dt\  Qc  j  dt 


c 


+  qV  ■ 


npi'l 

Qcj 


^  Pin-(n-V)v_^ 

~r  -)  —  U. 


(51) 


This  equation  coincides  with  Eq.  (3.18)  from  Ref.  9  if  in  the 
latter  equation  terms  proportional  to  V  -v  are  ignored.  Equa¬ 
tion  (3.18)  is  an  exact  transport  equation  for  pi  in  the  geo¬ 
metrical  acoustics  derived  from  Eqs.  (5)-(8).  Thus,  if  the 
terms  proportional  to  V-v  are  ignored,  Eqs.  (17)  and  (18) 
exactly  describe  the  amplitude  of  a  sound  wave  in  a  geo¬ 
metrical  acoustics  approximation,  and  correctly  account  for 
terms  of  any  order  in  vie.  Note  that  in  Ref.  9  starting  from 
the  transport  equation,  Eq.  (3.18),  a  law  of  acoustic  energy 
conservation  in  geometrical  acoustics  of  moving  media  is 
derived;  see  Eq.  (3.21)  from  that  reference.  Since  Eq.  (51) 
coincides  with  Eq.  (3.18),  the  same  law  [i.e.,  Eq.  (3.21)  from 
Ref.  9]  can  be  derived  from  Eq.  (51)  provided  that  the  terms 
proportional  to  V-v  are  ignored. 

G.  Discussion 

Thus,  by  comparing  a  set  of  Eqs.  (17)  and  (18)  with  the 
equations  for  p  which  are  widely  used  in  atmospheric  acous¬ 
tics,  we  determined  that  this  set  has  the  same  or  a  wider 
range  of  applicability  than  these  equations  for  p.  Note  that 
there  are  other  equations  for  p  known  in  the  literature  (see 
Refs.  9,  11,  17  and  references  therein):  Monin’s  equation, 
Pierce’s  equations,  equation  for  the  velocity  quasi-potential, 
the  Andreev-Rusakov-Blokhintzev  equation,  etc.  Most  of 
these  equations  have  narrower  ranges  of  applicability  than 
the  equations  presented  above  and  have  been  seldom  used 
for  calculations  of  p. 


V.  NUMERICAL  IMPLEMENTATION 


In  this  section,  we  describe  simple  algorithms  for  EDTD 
solutions  of  Eqs.  (17)  and  (18)  in  the  two  spatial  dimensions 
X  and  y.  Isolating  the  partial  derivatives  with  respect  to  time 
on  the  left  side  of  these  equations,  we  have 


dp 

Jt 


dy 


dw^  dw. 


dx 


dy 


|>H-. 

A 

•  p 


FIG.  1 .  Spatially  staggered  finite-difference  grid  used  for  the  calculations  in 
this  article. 


dt 


d  d\  d  d 

-b  —  +  bF^, 
dx 


(53) 


^  =  — ju.-ju,— +  i;,  — )w,, 


''dx'"ydy}y  [  ’‘dx  ^dyj  y 


dp 


(54) 


where  b  =  1/p  is  the  mass  buoyancy  and  k=  pc^  is  the  adia¬ 
batic  bulk  modulus.  In  Eqs.  (52)-(54),  the  subscripts  x  and  y 
indicate  components  along  the  corresponding  coordinate 
axes. 

The  primary  numerical  issues  pertinent  to  solving  these 
equations  in  a  moving  inhomogeneous  medium  are  summa¬ 
rized  and  addressed  in  Secs.  V  A-V  C.  Example  calculations 
are  provided  in  Secs.  V  D  and  VE. 


A.  Spatial  finite-difference  approximations 

The  spatial  finite-difference  (ED)  network  considered 
here  stores  the  pressure  and  particle  velocities  on  a  grid  that 
is  staggered  in  space,  as  shown  in  Eig.  1.  The  pressure  is 
stored  at  integer  node  positions,  namely  x  =  i  Ax  and  y 
=j  Ay,  where  i  and  j  are  integers  and  Ax  and  Ay  are  the 
grid  intervals  in  the  x-  and  y-directions.  The  x-components 
of  the  acoustic  velocity,  w^.,  are  staggered  (offset)  by  Ax/2 
in  the  x-direction.  The  y-components  of  the  acoustic  veloc¬ 
ity,  Wy,  are  staggered  by  Ay/2  in  the  y -direction.  This  stag¬ 
gered  grid  design  is  widely  used  for  wave  propagation  cal¬ 
culations  in  nonmoving  media. Here  we  furthermore 
store  and  at  the  nodes,  and  u,,  and  E,,  at  the  Wy 
nodes.  The  quantities  b,  k,  and  Q  are  stored  at  the  pressure 
nodes. 

Eor  simplicity,  we  consider  in  this  article  only  a  second- 
order  accurate,  spatially  centered  ED  scheme.  A  centered  so¬ 
lution  of  Eqs.  (52)-(54)  requires  an  evaluation  of  each  of  the 
terms  of  the  right-hand  sides  of  these  equations  at  the  grid 
nodes  where  the  field  variable  on  the  left-hand  side  is  stored. 
One  of  the  main  motivations  for  using  the  spatially  staggered 
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grid  is  that  it  conveniently  provides  compact,  centered  spatial 
differences  for  many  of  the  derivatives  in  Eqs.  (52)-(54). 
For  example,  dw^  /  dx  in  Eq.  (52)  is 

dw^{i  AxJ  Ay ,t)l dx—{w ^\_{i  +  \l2)Ax,j  Ay,f] 

—  w^\_{i—\l2)Ax,jAy,t~\}IAx  (55) 

and  dp! dy  in  Eq.  (54)  is 

dpl^i  Ax,(j+  l/2)Ay,f]/^y 

^{p[i  Ax,(j+l)Ay,t]-p[i  AxJ  Ay,f]}/Ay.  (56) 

The  derivatives  dpidx  and  dw^ldy  follow  similarly.  The 
body  source  terms  can  all  be  evaluated  directly,  since  they 
are  already  stored  at  the  grid  nodes  where  the  ED  approxi¬ 
mations  are  centered.  The  same  is  true  of  k,  which  is  stored 
at  the  pressure  grid  nodes  and  needed  in  Eq.  (52).  Regarding 
Eqs.  (53)  and  (54),  the  values  for  b  can  be  determined  at  the 
needed  locations  by  averaging  neighboring  grid  points. 

The  implementation  of  the  remaining  terms,  particular  to 
the  moving  medium,  is  somewhat  more  complicated.  For 
example,  the  derivatives  of  the  pressure  held  in  Eq.  (52), 
dpidx  and  dpidy,  cannot  be  centered  at  x  =  iAx  and  y 
=j  Ay  from  approximations  across  a  single  grid  interval. 
Centered  approximations  can  be  formed  across  two  grid  in¬ 
tervals,  however,  as  suggested  in  Ref.  2.  For  example, 

dp{i  AxJ  Ay,t)/dx 

—  {p[(;  -f  1  )AxJ  Ay,f]  —/?[(/—  1  )AxJ  Ay,f]}/2Ax. 

(57) 

Neighboring  grid  points  can  be  averaged  to  hnd  the  wind 
velocity  components  and  at  x  =  iAx  and  y=j  Ay, 
which  multiply  the  derivatives  dpidx  and  dpidy,  respec¬ 
tively,  in  Eq.  (52).  Similarly,  the  spatial  derivatives  of  the 
particle  velocities  in  Eqs.  (53)  and  (54)  can  be  approximated 
over  two  grid  intervals.  In  Eq.  (53),  the  quantities  Wy  and  Vy 
(multiplying  the  derivatives  dv^ldy  and  dw^ldy,  respec¬ 
tively)  are  needed  at  the  grid  point  x  =  (; -f  1/2) Ax  and  y 
=j  Ay.  Referring  to  Fig.  1,  a  reasonable  way  to  obtain  these 
quantities  would  be  to  average  the  four  closest  grid  nodes; 

^,,[(1  +  112) Ax, j  Ay,f] 

^  4  +  1  1/2) Ay, f] 

-f-Wy[i  Ax,{j+  l/2)Ay,f] 

+  w^,[(/+  l)Ax,(;-  l/2)Ay,f] 

+  Wj,[i  Ax,(;-l/2)Ay,f]},  (58) 

and  likewise  for  Vy.  The  quantities  and  v^,  multiplying 
the  derivatives  dvyidx  and  dWyldx  in  Eq.  (54),  can  be  ob¬ 
tained  similarly. 


B.  Advancing  the  solution  in  time 

Let  us  dehne  the  functions  fp,  f^,  and  fy  as  the  right- 
hand  sides  of  Eqs.  (52),  (53),  and  (54),  respectively.  For 
example,  we  write 


dpj  AxJ  Ay,t) 

Jt 

=fp[i  AxJ  Ay,p(f),W;,(f),Wj,(f),s(f)],  (59) 

where  p(f),  ’wjt),  and  Wj,(f)  are  matrices  containing  the 
pressures  and  acoustic  velocities  at  all  available  grid  nodes. 
For  convenience,  s(f)  is  used  here  as  short  hand  for  the  com¬ 
bined  source  and  medium  properties  (b,  k,  v^,  Vy,  Q,  F^, 
and  Fy)  at  all  available  grid  nodes.  (Note  that 
fp[iAx,jAy,p(t),wJt),-Wy{t),s(t)]  in  actuality  depends 
only  on  the  helds  at  a  small  number  of  neighboring  grid 
points  of  {i  Ax,j  Ay)  when  second-order  spatial  differencing 
is  used.  The  notation  here  is  general  enough,  though,  to  ac¬ 
commodate  spatial  differencing  of  an  arbitrarily  high  order.) 

For  a  nonmoving  medium,  the  solution  is  typically  ad¬ 
vanced  in  time  using  a  staggered  temporal  grid,  in  which  the 
pressures  are  stored  at  the  integer  time  steps  t  =  l  At  and  the 
particle  velocities  at  the  half-integer  time  steps  t  =  {l 
+  l/2)Af.^^“^°  The  acoustic  velocities  and  pressures  are  up¬ 
dated  in  an  alternating  “leap-frog”  fashion,  with  the  helds 
from  the  previous  time  step  being  overwritten  in  place.  Con¬ 
sidering  the  moving  media  equations,  approximation  of  the 
time  derivative  in  Eq.  (59)  with  a  hnite  difference  centered 
on  f  =  (/-I- 1/2) Af  (that  is,  <?/?[/ Ax,y  Ay, (/-f  l/2)Af]/(?f 
—  {p[;  Ax,i  Ay,(/-I-  l)Af]— p[/  Ax,]  Ay, I  Af]}/Af)  results 
in  the  following  equation  for  updating  the  pressure  held: 

p[i  Ax,j  Ay,{l+  1 ) Af] 

=p[i  AxJ  Ay, I  Af]-l-Af/^,[/  AxJ  Ay,p[(/-I- 1/2) Af], 

wi(/+  l/2)Af],W3,[(/-f  l/2)Af],s[(/-l-  l/2)Af]]. 

(60) 

Note  that  this  equation  requires  the  pressure  held  at  the  half¬ 
integer  time  steps,  i.e.,  f  =  (/-!- 1/2) Af.  In  the  staggered  leap¬ 
frog  scheme,  however,  the  pressure  is  unavailable  at  the  half¬ 
integer  time  steps.  A  similar  centered  approximation  for  the 
acoustic  velocities  indicates  that  they  are  needed  on  the  in¬ 
teger  time  steps  in  order  to  advance  the  solution,  which  is 
again  problematic.  If  one  attempts  to  address  this  problem  by 
linearly  interpolating  between  adjacent  time  steps  (i.e.,  by 
setting  p[(/-l- l/2)Af]  — {p[/ Af]-l-p[(/-l- l)Af]}/2  in  Eq. 
(60)),  explicit  updating  equations  (a  solution  of  Eq.  (60)  for 
/?[/ Ax,y  Ay,(/-I- 1 ) Af]  that  does  not  require  the  pressure 
held  at  nearby  grid  points  at  the  time  step  t  =  {l+  l)Af)  can¬ 
not  be  obtained.  Hence  the  customary  staggered  leap-frog 
approach  does  not  lead  to  an  explicit  updating  scheme  for  the 
acoustic  helds  in  a  moving  medium.  The  staggered  leap-frog 
scheme  can  be  rigorously  implemented  only  when  the  terms 
particular  to  the  moving  medium  (those  involving  and  Vy) 
are  removed  from  Eqs.  (52)-(54). 

A  possible  work-around  would  be  to  use  the  pressure 
held  p(/  Af)  in  place  of  p[(/  -f  1/2) Af]  when  evaluating  fp , 
and  Wjj[(/ — 1/2) Af]  and  w,,[(/— 1/2) Af]  in  place  of 
w^(/  Af)  and  Wj,(/  Af)  when  evaluating /j.  and  fy  .  This  non¬ 
rigorous  procedure  uses  the  Euler  (forward  difference) 
method  to  evaluate  the  moving-media  terms  while  maintain¬ 
ing  the  leap-frog  approach  for  the  remaining  terms.  From  a 
programming  standpoint,  the  algorithm  proceeds  in  essen- 
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tially  the  same  manner  as  the  staggered  leap-frog  method  for 
a  nonmoving  medium.  The  calculations  in  Ref  2  appear  to 
use  such  a  procedure.  But  the  stability  and  accuracy  of  this 
algorithm  are  unclear.  An  alternative  is  provided  in  Ref.  4, 
which  uses  a  perturbative  solution  based  on  the  assumption 
that  the  flow  velocity  is  small. 

Here  we  would  like  to  develop  a  general  technique  that 
is  applicable  to  high  Mach  numbers.  The  simplest  way  to 
accomplish  this  is  to  abandon  the  staggered  temporal  grid 
and  form  centered  finite  differences  over  two  time  steps. 
The  pressure  updating  equation,  based  on  the  approximation 
dp(i  Ax,  j  Ay,  I  At)/ dt—{p[i  Ax,  /  Ay,(/+ l)Af] 
-p[i  Ax,  jAy,{l-l)At]}/2At,  is 

p[i  Ax,j  Ay, {I  +  1  )Af] 

=p[i  Ax,j  Ay,{[-l)At] 

+  2  Atfp[i  Ax,j  Ay,pil  At),wJJ  At), 

Wyil  At),s{l  At)].  (61) 

Similarly,  we  derive 

+  1/2) Ax, j  Ay,{l+  l)Af] 

=  w^[{i+ 1/2) Ax, j  Ay,{l-l)At]  +  2Atfj,[{i 

+  1/2) Ax, j  Ay,p(l  At),w^{l  At),Wy{l  At),s{l  Af)], 

(62) 

Ax,(j+  1/2) Ay, {I  +  l)Af] 

=  Wy[i  Ax,{j+l/2)Ay,(l-l)At] 

+  2  Atfy[i  Ax,{j+l/2)Ay,p(l  Af),w^ 

(I  At),w,{l  At),s{l  At)].  (63) 

Somewhat  confusingly,  this  general  temporal  updating 
scheme  has  also  been  called  the  “leap-frog”  scheme  in  the 
literature,^*  since  it  involves  alternately  overwriting  the 
wavefield  variables  at  even  and  odd  integer  time  steps  based 
on  calculations  with  the  fields  at  the  intervening  time  step. 
We  call  this  scheme  here  the  nonstaggered  leap-frog.  The 
primary  disadvantage,  in  comparison  to  the  staggered  leap¬ 
frog  scheme,  is  that  the  fields  must  be  stored  over  two  time 
steps,  rather  than  just  one.  Additionally,  the  numerical  dis¬ 
persion  and  instability  characteristics  are  inferior  to  those  of 
the  conventional  staggered  scheme  due  to  the  advancement 
of  the  wavefield  variables  over  two  time  steps  instead  of  one. 
On  the  other  hand,  the  nonstaggered  leap-frog  does  provide  a 
simple  and  rigorous  centered  finite-difference  scheme  that  is 
not  specialized  to  low  Mach  number  flows.  Other  common 
numerical  integration  methods,  such  as  the  Runge-Kutta 
family,  can  also  be  readily  applied  to  the  nonstaggered-in- 
time  grid.  Some  of  the  calculations  following  later  in  this 
section  use  a  fourth-order  Runge-Kutta  method,  which  is 
described  in  Ref.  32  and  many  other  texts.  We  have  also 
developed  a  staggered-in-time  method  that  is  valid  for  high 
Mach  numbers  but  requires  the  fields  to  be  stored  over  two 
time  levels.  This  method  was  briefly  discussed  in  Ref.  6. 

Note  that  our  present  numerical  modeling  efforts  are  di¬ 
rected  toward  demonstrating  the  applicability  and  feasibility 


of  FDTD  techniques  for  simulating  sound  propagation  in  a 
moving  atmosphere.  We  have  not  undertaken  a  comprehen¬ 
sive  comparative  analysis  of  the  many  alternative  numerical 
strategies  available  for  the  solution  of  Eqs.  (17)  and  (18). 
However,  several  of  these  approaches  (including  the  pseu- 
dospectral  method,  higher-order  spatial  and/or  temporal 
finite-difference  operators,  and  the  dispersion  relation  pre¬ 
serving  (DRP)  technique)  yield  accurate  simulations  of 
sound  propagation  with  fewer  grid  intervals  per  wavelength 
compared  with  our  numerical  examples.  In  particular,  the 
DRP  method,  involving  optimized  numerical  values  of  the 
finite-difference  operator  coefficients  (e.g..  Ref.  18),  can  be 
readily  introduced  into  our  FDTD  algorithmic  framework. 


C.  Dependence  of  grid  increments  on  Mach  number 

For  numerical  stability  of  the  2-D  FDTD  calculation,  the 
time  step  Af  and  grid  spacing  Ar  must  be  chosen  to  satisfy 
the  Courant  condition,  C<  l/v5  (e.g.,  see  Ref.  33),  where  the 
Courant  number  is  defined  as 

u  At 


Here,  u  is  the  speed  at  which  the  sound  energy  propagates. 
[For  a  nonuniform  grid,  Ar=  l/V(Ax)^^-l-(Ay)^^.]  Since 
the  grid  spacing  must  generally  be  a  small  fraction  of  a 
wavelength  for  good  numerical  accuracy,  the  Courant  condi¬ 
tion  in  practice  imposes  a  limitation  on  the  maximum  time 
step  possible  for  stable  calculations.  An  even  smaller  time 
step  may  be  necessary  for  good  accuracy,  however. 

Let  us  consider  the  implications  of  the  Courant  condi¬ 
tion  for  propagation  in  a  uniform  flow.  In  this  case,  u  is 
determined  by  a  combination  of  the  sound  speed  and  wind 
velocity.  In  the  downwind  direction,  we  have  m  =  m  +  =  c 
-f  u .  In  the  upwind  direction,  «  =  «  _  =  c  —  u .  The  wave¬ 
lengths  in  these  two  directions  are  \^  =  {c  +  v)/ f  and  X._ 
=  {c  —  v)/f,  respectively,  where  /  is  the  frequency.  Since  the 
wavelength  is  shortest  in  the  upwind  direction,  the  value  of 
X._  dictates  the  grid  spacing.  We  set 


where  N  is  the  number  of  grid  points  per  wavelength  in  the 
upwind  direction,  M  =  v/c  is  the  Mach  number,  and  k  =  c/f 
is  the  wavelength  for  the  medium  at  rest.  If  N  is  to  be  fixed 
at  a  constant  value,  a  finer  grid  is  required  as  M  increases. 
Regarding  the  time  step,  the  Courant  condition  implies 

Af<^.  (66) 

This  condition  is  most  difficult  to  meet  when  u  is  largest, 
which  is  the  case  in  the  downwind  direction.  Therefore  we 
must  use  in  the  preceding  inequality  if  we  are  to  have 
accurate  results  throughout  the  domain;  specifically,  we  must 
set 


At< 


nu7 


1  1-M 

j/f  r+M' 


(67) 
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FIG.  2.  The  geometry  of  the  problem. 


Therefore  the  time  step  must  also  be  shortened  as  M  in¬ 
creases.  For  example,  the  time  step  at  M=  1/3  must  be  1/2 
the  value  necessary  at  M  =  0.  At  M  =  2/3,  the  time  step  must 
be  1/5  the  value  at  M  =  0.  The  reduction  of  the  required  time 
step  and  grid  spacing  combine  to  make  calculations  at  large 
Mach  numbers  computationally  expensive. 


D.  Example  calculations 

In  this  subsection,  we  use  the  developed  algorithm  for 
FDTD  solutions  of  Eqs.  (52)-(54)  to  compute  the  sound 
field  p  in  a  2-D  homogeneous  uniformly  moving  medium. 
The  geometry  of  the  problem  is  shown  in  Fig.  2.  A  point 
monochromatic  source  is  located  at  the  origin  of  the  Carte¬ 
sian  coordinate  system  x,y.  The  medium  velocity  v  is  paral¬ 
lel  to  the  x-axis.  We  will  first  obtain  an  analytical  formula  for 
p  for  this  geometry. 

In  a  homogeneous  uniformly  moving  medium,  c,  Q,  and 
V  are  constant  so  that  V-v=0  and  VP  =  Q.  Therefore,  Eqs. 
(17)  and  (18)  describe  sound  propagation  exactly  for  this 
case  and  are  valid  for  an  arbitrary  value  of  the  Mach  number 
M.  They  can  be  written  as 


— -FvV  \p  +  ■yi=  Qc^Q, 

(68) 

d  \  Vp 

—  -F  vV  w-l - =  0. 

dt  j  Q 

(69) 

Here,  p  and  w  are  functions  of  the  coordinates  x,  y  and  time 
t,  V  =  {d/dx,d/dy),  and  the  function  Q  is  given  by 

2iA 

Q=—e-‘'^‘S(x)S{y),  (70) 

^  0) 

where  S  is  the  delta  function  and  the  factor  A  characterizes 
the  source  amplitude.  In  Eqs.  (68)  and  (69),  for  simplicity,  it 
is  assumed  that  F  =  0. 

Assuming  that  v<c,  the  following  solution  of  Eqs.  (68) 
and  (69)  is  obtained  in  the  Appendix: 

p(r,a,M) 


iA 


2(1-M^) 


2^3/2 


iM  cos  a 
Vl  —M^  sin^  a 


(^) 


Xexp 


ikMr  cos  a 
1  -M^ 


(71) 


Here,  k=u)/c,  and  H[^'^  are  the  Hankel  functions,  ^ 

=  krVl  —M^  sin"'  and  r  and  a  are  the  polar  coor¬ 

dinates  shown  in  Eig.  2.  Eor  krS>  1,  the  Hankel  functions  can 
be  approximated  by  their  asymptotics.  This  results  in  the 
desired  formula  for  the  sound  pressure: 

A(  Vl  —M^  sin^  a  —  M  cos  a) 

pir  OL  A/)  —  — ^ - 

V27r/tr(l -M2)(1 sin^  a)3/4 


Xexp 


sin^  a  —  M  cos  a)kr  iir 
l-M^ 


(72) 


Note  that  a  sound  field  due  to  a  point  monochromatic  source 
in  a  2-D  homogeneous  uniformly  moving  medium  was  also 
studied  in  Ref.  18  by  a  different  approach.  The  phase  factor 
obtained  in  that  reference  is  essentially  the  same  as  that  in 
Eq.  (72).  Only  a  general  expression  for  the  amplitude  factor 
was  presented  in  Ref.  18  which  does  not  allow  a  detailed 
comparison  with  the  amplitude  factor  in  Eq.  (72). 

Let  us  now  consider  the  EDTD  calculations  of  the  sound 
field  for  the  geometry  in  Eig.  2.  In  these  calculations,  the 
source  consists  of  a  finite-duration  harmonic  signal  with  a 
cosine  taper  function  applied  at  the  beginning  and  the  end. 
The  tapering  alleviates  numerical  dispersion  of  high  frequen¬ 
cies,  which  becomes  evident  when  there  is  an  abrupt  change 
in  the  source  emission.  The  tapered  source  equation  is 

(1/2)[1  — cos(T7f/7’i)]cos(27r/-l-  <p), 

0«f<7’i, 

cos{2 TTf  +  (})) ,  T’lSSfssT’— ^2, 

(1/2)[1  -l-cos(7r(f  — 7’)/7’2)]cos(27r/-l-  (/>), 
T-T2<t^T, 

0,  otherwise. 


G(0  = 


(73) 


Here,  </>  is  the  source  phase,  T j  is  the  duration  of  the  initia¬ 
tion  taper,  and  7’2  is  the  duration  of  the  termination  taper.  All 
calculations  in  this  paper  use  tapering  over  an  interval  of  3 
periods  in  the  harmonic  wave  ( Tj  =  7’2  =  3// )  and  a  total 
signal  duration  of  10  periods  {T=  10//). 

Eigure  3  shows  the  pressure  field  for  a  100  Hz  source  in 
a  uniform  Mach  0.3  flow.  The  field  is  shown  at  0.11  s,  or 
0.01  s  after  the  source  has  been  turned  off.  The  distance 
between  wave  fronts  is  smaller  upwind  than  downwind.  The 
calculations  use  the  fourth-order  Runge-Kutta  method  with 
a  staggered  spatial  grid  and  a  nonstaggered  temporal  grid. 
The  spatial  domain  is  100  m  by  100  m,  with  800  grid  points 
in  each  direction.  This  results  in  approximately  19  grid 
points  per  wavelength  in  the  upwind  direction.  The  time  step 
was  set  to  0.145  ms,  which  implies  a  Courant  number  of  0.40 
in  a  nonmoving  medium  but  0.52  in  the  downwind  direction 
of  the  M  =  0.3  flow.  Using  the  run  shown  in  Eig.  3,  the  azi¬ 
muthal  dependence  of  the  normalized  sound  pressure  magni¬ 
tude  \p{r,a,M)lp{rfl,0)\  for  values  of  kr  ranging  from  1  to 
100  was  compared  to  the  theoretical  far-fleld  result  calcu¬ 
lated  from  Eq.  (72).  Excellent  agreement  was  found  between 
theoretical  predictions  and  EDTD  simulations  for  krSlO. 

The  azimuthal  dependence  of  |/2(r,Q',M)//2(r,0,0)|  at 
kr  =  20  is  compared  for  several  numerical  methods  in  Eig.  4. 
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FIG.  3.  Wavefronts  of  the  sound  pressure  due  to  a  point 
source  located  at  the  point  x  =  0  and  v  =  0  for  M 
=  0.3.  The  medium  velocity  is  in  the  direction  of  the 
x-axis. 


The  methods  include  the  staggered  (with  forward- 
differencing  of  the  moving  medium  terms  mentioned  in  Sec. 
V  B)  and  nonstaggered  leap-frog  approaches  and  the  fourth- 
order  Runge-Kutta.  The  time  step  for  the  leap-frog  methods 
was  0.036  ms  (1/4  that  used  for  the  Runge-Kutta),  so  that 
the  computational  times  of  all  calculations  are  roughly  equal. 
The  Runge-Kutta  and  nonstaggered  leap-frog  provide 
graphically  indistinguishable  results.  The  staggered  leap¬ 
frog,  however,  systematically  underpredicts  the  amplitude  in 
the  downwind  direction  and  overpredicts  in  the  upwind  di¬ 
rection.  The  actual  sound  pressure  signals  at  f  =  0. 1 1  s,  cal¬ 
culated  from  the  staggered  and  nonstaggered  leap-frog  ap¬ 
proaches,  are  overlaid  in  Fig.  5.  In  the  downwind  direction, 
the  staggered  leap-frog  method  provides  a  smooth  prediction 
at  distances  greater  than  about  22  m.  The  noisy  appearance  at 
shorter  distances  is  due  to  numerical  instability,  which  was 
clear  from  the  rapid  temporal  growth  of  these  features  we 


observed  as  the  calculation  progressed.  We  conclude  that  the 
staggered  leap-frog  approach,  when  applied  to  a  moving  me¬ 
dium,  is  less  accurate  and  more  prone  to  numerical  instabil¬ 
ity.  This  is  likely  due  to  the  nonsymmetric  temporal  finite 
difference  approximations  for  the  moving  medium  terms. 

Figure  6  shows  the  azimuthal  dependence  of 
\p{r,a,M)lp{r,0,0)\  for  M  =  0,  0.3,  and  0.6.  All  FDTD  cal¬ 
culations  for  this  figure  use  the  fourth-order  Runge-Kutta 
method.  Two  calculated  curves  are  shown;  one  for  a  low- 
resolution  run  with  800X800  grid  points  and  a  time  step  of 
0.145  ms,  and  the  other  for  a  high-resolution  run  with  1600 
X  1600  grid  points  and  a  time  step  of  0.0362  ms.  For  M 
=  0.3,  both  grid  resolutions  yield  nearly  exact  agreement 
with  Eq.  (72).  At  M  =  0.6,  the  low-resolution  run  has  11 
spatial  grid  nodes  per  wavelength  in  the  upwind  direction 
and  a  downwind  Courant  number  of  0.64.  The  high- 
resolution  grid  has  22  spatial  grid  nodes  per  wavelength  in 


FIG.  4.  Normalized  sound  pressure  amplitude 
\p{r  , a, M)/p(r, 0,0)1  versus  the  azimuthal  angle  a  for 
M  =  0.3  and  kr  =  2Q.  The  staggered  and  nonstaggered 
leap-frog  methods  and  the  fourth-order  Runge-Kutta 
are  compared  to  the  theoretical  solution.  The  nonstag¬ 
gered  leap-frog  and  Runge-Kutta  methods  are  graphi¬ 
cally  indistinguishable. 
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Downwind 


Upwind 


FIG.  5.  Sound  pressure  traces  for  (a)  downwind  and  (b) 
upwind  propagation.  Calculations  from  the  staggered 
and  nonstaggered  leap-frog  methods  are  shown  (dashed 
and  solid  lines,  respectively). 


the  upwind  direction  and  a  downwind  Courant  number  of 
0.32.  Agreement  with  theory  at  M  =  0.6  is  very  good  for  the 
high-resolution  run.  The  low-resolution  run  substantially  un¬ 
derpredicts  the  upwind  amplitude. 

Finally  note  that  it  follows  from  Figs.  4  and  6  that  the 
sound  pressure  is  largest  for  a'=180°,  i.e.,  in  the  upwind 
direction.  This  dependence  is  also  evident  upon  close  exami¬ 
nation  in  Fig.  3. 

E.  Comparison  of  FDTD  and  FFP  calculations 

The  computational  examples  so  far  in  this  paper  have 
been  for  uniform  flows.  However,  the  numerical  methods 
and  equations  upon  which  they  are  based  apply  to  nonuni¬ 
form  flows  as  well.  In  this  section,  we  consider  an  example 
calculation  for  a  flow  with  constant  shear.  The  point  source 
and  receiver  are  both  located  at  a  height  of  20  m  and  the 
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frequency  is  100  Hz.  The  computational  domain  is  200  m  by 
100  m  and  has  600  by  300  grid  points.  The  time  step  is 
7.73X10^^  s  and  the  fourth-order  Runge-Kutta  method  is 
used.  A  rigid  boundary  condition  is  applied  at  the  ground 
surface  (y  =  0  m).  An  absorbing  layer  in  the  upper  one-fifth 
of  the  simulation  domain  removes  unwanted  numerical  re¬ 
flections.  (The  implementation  of  the  rigid  ground  boundary 
condition  and  the  absorbing  layer  is  described  in  Ref.  34. 
Realistic  ground  boundary  conditions  in  a  FDTD  simulation 
of  sound  propagation  in  the  atmosphere  are  considered  in 
Ref.  35.) 

Calculated  transmission  loss  (sound  level  relative  to  free 
space  at  1  m  from  the  source)  results  are  shown  in  Figs.  7(a) 
and  7(b).  The  first  of  these  figures  is  for  a  zero-wind  condi¬ 
tion  and  the  second  is  for  a  horizontal  (x-direction)  wind 
speed  of  v(y)  =  fJiy,  where  the  gradient  /x  is  1  s^*.  For  Fig. 
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FIG.  6.  Normalized  sound  pressure  amplitude 
\pir  ,a,M)//7(r,0,0)|  versus  the  azimuthal  angle  a  for 
M  =  0,  0.3,  and  0.6.  The  fourth-order  Runge-Kutta 
method  was  used.  The  calculation  with  800X800  grid 
points  had  a  spatial  resolution  of  0.125  m  and  time  step 
0.145  ms,  whereas  the  1600X1600  calculation  had  a 
spatial  resolution  of  0.0625  m  and  time  step  0.0362  ms. 


7(a),  the  FDTD  results  are  compared  with  both  the  exact 
solution  for  a  point  source  above  the  rigid  boundary  and 
calculations  from  the  FFP  developed  in  Ref.  36.  The  FDTD 
results  are  nearly  indistinguishable  from  the  exact  solution. 
The  FFP  is  also  in  good  agreement,  although  there  is  some 
systematic  underprediction  of  the  interference  minima,  par¬ 
ticularly  so  near  the  source.  This  is  likely  due  to  the  far-held 
approximation  inherent  to  the  FFP.  For  the  case  with  con¬ 
stant  shear.  Fig.  7(b),  the  interference  pattern  is  shifted.  The 
FDTD  and  FFP  continue  to  show  very  similar  small  discrep¬ 
ancies  near  the  source.  On  the  basis  of  the  results  shown  in 
Fig.  7(a),  it  is  highly  likely  that  the  FDTD  is  more  accurate. 
The  FDTD  calculations  required  about  100  times  as  long  to 
complete  as  the  FFP  on  a  single-processor  computer.  As 
would  be  expected,  the  FFP  is  more  efficient  for  calculations 
at  a  limited  number  of  frequencies  in  a  horizontally  stratihed 
medium. 

VI.  CONCLUSIONS 

In  the  present  paper,  we  have  considered  starting  equa¬ 
tions  for  FDTD  simulations  of  sound  propagation  in  a  mov¬ 
ing  inhomogeneous  atmosphere.  FDTD  techniques  can  pro¬ 
vide  a  very  accurate  description  of  sound  propagation  in 
complex  environments. 

A  most  general  description  of  sound  propagation  in  a 
moving  inhomogeneous  medium  is  based  on  the  complete 
set  of  linearized  equations  of  fluid  dynamics,  Eqs.  (5)-(8). 
However,  this  set  is  too  involved  to  be  effectively  employed 
in  FDTD  simulations  of  outdoor  sound  propagation.  In  this 
paper,  the  linearized  equations  of  fluid  dynamics  were  re¬ 
duced  to  two  simpler  sets  of  equations  which  can  be  used  as 
starting  equations  for  FDTD  simulations. 

The  first  set  of  equations  contains  three  coupled  equa¬ 
tions,  Eqs.  (5),  (6),  and  (13),  for  the  sound  pressure  p,  acous¬ 
tic  velocity  w,  and  acoustic  density  rj.  This  set  is  an  exact 
consequence  of  the  linearized  equations  of  fluid  dynamics, 
Eqs.  (5)-(8).  To  solve  the  first  set  of  equations,  one  needs  to 
know  the  following  ambient  quantities:  the  adiabatic  sound 


speed  c,  density  Q,  medium  velocity  v,  pressure  P,  and  the 
parameters  a,  /8,  and  h .  The  atmosphere  can  be  modeled  as 
an  ideal  gas  to  a  very  good  accuracy.  In  this  case,  the  first  set 
of  equations  simplifies  and  is  given  by  Eqs.  (5),  (6),  and  (15). 
Now  it  contains  the  following  ambient  quantities:  c,  Q,  v, 
and  P. 

The  second  set  of  starting  equations  for  EDTD  simula¬ 
tions  contains  two  coupled  equations  for  the  sound  pressure 
p  and  acoustic  velocity  w,  Eqs.  (17)  and  (18).  In  order  to 
solve  this  set  one  needs  to  know  a  fewer  number  of  the 
ambient  quantities:  c,  Q,  and  v.  Note  that  namely  these  am¬ 
bient  quantities  appeared  in  most  of  equations  for  the  sound 
pressure  p  which  have  been  previously  used  for  analytical 
and  numerical  studies  of  outdoor  sound  propagation.  The 
second  set  was  derived  from  Eqs.  (5)-(8)  assuming  that 
terms  proportional  to  the  divergence  of  the  medium  velocity 
and  the  gradient  of  the  ambient  pressure  can  be  ignored. 
Both  these  assumptions  are  reasonable  in  atmospheric  acous¬ 
tics.  To  better  understand  the  range  of  applicability  of  the 
second  set,  it  was  compared  with  equations  for  the  sound 
pressure  p  which  have  been  most  often  used  for  analytical 
and  numerical  studies  of  sound  propagation  in  a  moving  in¬ 
homogeneous  medium.  It  was  shown  that  the  second  set  has 
the  same  or  wider  range  of  applicability  than  these  equations 
for  p.  Thus,  a  relatively  simple  set  of  Eqs.  (17)  and  (18), 
which  is  however  rather  general,  seems  very  attractive  as 
starting  equations  for  EDTD  simulations. 

The  numerical  algorithms  for  EDTD  solutions  of  the 
second  set  of  equations  were  developed  for  the  case  of  a  2-D 
inhomogeneous  moving  medium.  It  was  shown  that  the 
staggered-in-time  grid  approach  commonly  applied  to  non¬ 
moving  media  cannot  be  applied  directly  for  the  moving 
case.  However,  fairly  simple  alternatives  based  on 
nonstaggered-in-time  grids  are  available.  We  used  the  result¬ 
ing  algorithms  to  calculate  the  sound  pressure  due  to  a  point 
source  in  a  homogeneous  uniformly  moving  medium.  The 
results  obtained  were  found  in  excellent  agreement  with  ana¬ 
lytical  predictions  even  for  a  Mach  number  as  high  as  0.6. 
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FIG.  7.  Comparisons  between  the  transmission  loss  calculated  with  different 
methods,  (a)  Homogeneous  atmosphere  without  wind,  (b)  Atmosphere  with 
linearly  increasing  wind  velocity. 


Furthermore,  using  the  algorithm  developed,  we  calculated 
the  sound  field  due  to  a  point  source  in  a  stratified  moving 
atmosphere.  The  results  obtained  are  in  a  good  agreement 
with  the  FFP  solution. 

Finally  note  that  Eqs.  (17)  and  (18)  have  already  been 
used  as  starting  equations  in  FDTD  simulations  of  sound 
propagation  in  3-D  moving  media  with  realistic  velocity 
fields.  The  results  obtained  were  published  in  proceedings  of 
conferences. These  realistic  velocity  fields  include  the  fol¬ 
lowing:  kinematic  turbulence  generated  by  quasi-wavelets,^’® 
3-D  stratihed  moving  atmosphere,®  and  atmospheric  turbu¬ 
lence  generated  by  large-eddy  simulation.^  In  Ref.  8,  FDTD 
simulations  were  used  to  numerically  study  infrasound 
propagation  in  a  moving  atmosphere  over  distances  of  sev¬ 
eral  hundred  km.  The  largest  run  to  date  incorporated  over 
1.5  billion  nodes  and  took  about  100  hours  on  500  Compaq 
EV6  parallel  processors.* 
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APPENDIX:  SOUND  FIELD  DUE  TO  A  POINT 
MONOCHROMATIC  SOURCE  IN  A  HOMOGENEOUS 
UNIFORMLY  MOVING  MEDIUM 


In  this  appendix,  we  derive  a  formula  for  the  sound  pres¬ 
sure  due  to  a  point  monochromatic  source  located  in  a  2-D 
homogeneous  uniformly  moving  medium  (see  Fig.  2). 

For  this  geometry,  Eqs.  (68)  and  (69)  can  be  reduced  to 
a  single  equation  for  the  sound  pressure: 

(7+,.vj  p-c=vV=ec=(4+,.v|e.  (Ai) 

Here,  the  source  function  Q  is  given  by  Eq.  (70)  and  contains 
the  time  factor  exp(— /wf).  In  what  follows,  this  time  factor  is 
omitted.  Furthermore,  taking  into  account  that  the  medium 
velocity  is  parallel  to  the  x-axis,  Eq.  (Al)  can  be  written  as 

I  ,  d  .  d^\ 

-T-^  +  k^  +  likM  - - M^—^\p{x,y) 

\dx~  dy^  dx  dx^  I 

2iA  I  d  \ 

= -  /w-u—  5(x)(5(y).  (A2) 

w  \  dx  I 

Let 


2iA  /  d 

p{x,y)=  iut-v  —  \^{x,y). 


(A3) 


Substituting  this  formula  into  Eq.  (A2),  we  obtain  the  follow¬ 
ing  equation  for  the  function  <I>(x,y): 

<I)(x,y)  =  (5(x)5(y). 

(A4) 

In  this  equation,  let  us  make  the  following  transformations: 
x=  Vl-M^Y,  k=^\-M^K, 


^2 

d^ 

dY] 

- T  H“ 

-ik  +  M  —  \ 

dx^ 

dx] 

d)(x,y)  =  exp(-//:MA)^(A,y). 


(A5) 


As  a  result,  we  obtain  the  following  equation  for  the  fuinc- 
tion  ^(A,y): 


d^ 


d^ 


dX^  dy^ 


■+K- 


'i'iX,y)  = 


1 


Vl-M- 


:<5(Y)5(y).  (A6) 


A  solution  of  this  equation  is  well  known: 


'P(A,y) 


i 


H\l\K^X^  +  y^). 


(A7) 


Using  this  expression  for  "'P  and  Eqs.  (A3)  and  (A5),  we 
obtain  a  desired  formula  for  the  sound  pressure  of  a  point 
monochromatic  source  in  a  2-D  homogeneous  uniformly 
moving  medium: 


p(x,y)  = 


iA 


2(1-M^) 


2^3/2 


iMkx 


/  ixkM  \ 

(A8) 
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Here,  —M^)  +y^.  In  polar  coordi¬ 

nates,  Eq.  (A8)  becomes  Eq.  (71). 
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